/** @file fftw_impl.h
 * @brief Implementation of fftw-enabled FFT engines
 *
 */
#ifndef QMCPLUSPLUS_FFTW_IMPL_H
#define QMCPLUSPLUS_FFTW_IMPL_H
#include <fftw3.h>

//namespace fft_v0 {

/** base class for FFT engines using fftw library
 */
struct fftw_base
{
  typedef fftw_plan fft_plan_type;
  fft_plan_type forward_plan;
  fft_plan_type backward_plan;

  ///default constructor
  fftw_base() :forward_plan(0), backward_plan(0)
  { }

  ///virtual destructor to clean up the plans
  virtual ~fftw_base()
  {
#pragma omp critical
    {
      clean();
    }
  }

  void clean()
  {
    if(forward_plan) {fftw_destroy_plan(forward_plan); forward_plan=0;}
    if(backward_plan) {fftw_destroy_plan(backward_plan); backward_plan=0;}
  }

  template<typename T1>
  inline fftw_complex* magle(T1* in)
  {
    return reinterpret_cast<fftw_complex*>(in);
  }
};


/** specialization for single-fft only with FFTW on complex
 */
template<typename T>
struct FFT_DFT<T,std::complex<T>,false,FFTW_ENG> 
: public fft_base, public fftw_base
{
  typedef std::complex<T> value_type;

  FFT_DFT(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    init(dims,in,1,uflag);
  }

  void init(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    if(dims != fft_size)
    {
      fft_base::resize(dims,1);
#pragma omp critical
      {
        fftw_base::clean();
        forward_plan = fftw_plan_dft_1d(dims,magle(in),magle(in),FFTW_FORWARD,uflag);
        backward_plan = fftw_plan_dft_1d(dims,magle(in),magle(in),FFTW_BACKWARD,uflag);
      }
    }
  }

  inline void transform_forward_atom(value_type* in)
  {
    fftw_execute_dft(forward_plan,magle(in),magle(in));
  }

  inline void transform_backward_atom(value_type* in)
  {
    fftw_execute_dft(backward_plan,magle(in),magle(in));
  }

  inline void transform_forward(value_type* in, int howmany=1)
  {
    value_type* offset=in;
    for(int i=0; i<howmany; ++i,offset+=fft_size) 
      transform_forward_atom(offset);
  }

  inline void transform_backward(value_type* in, int howmany=1)
  {
    value_type* offset=in;
    for(int i=0; i<howmany; ++i,offset+=fft_size) 
      transform_backward_atom(offset);
  }
};

/** specialization for single-fft only with FFTW on double
 */
template<typename T>
struct FFT_DFT<T,T,false,FFTW_ENG> 
: public fft_base, public fftw_base
{
  typedef T value_type;
  int fft_mem_offset;

  FFT_DFT(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    init(dims,in,1,uflag);
  }

  void init(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    if(dims != fft_size)
    {
      fft_mem_offset=dims+2;
      fft_base::resize(dims,1);
#pragma omp critical
      {
        fftw_base::clean();
        forward_plan = fftw_plan_dft_r2c_1d(dims,in,magle(in),uflag);
        backward_plan = fftw_plan_dft_c2r_1d(dims,magle(in),in,uflag);
      }
    }
  }

  inline void transform_forward_atom(value_type* in)
  {
    fftw_execute_dft_r2c(forward_plan,in,magle(in));
  }

  inline void transform_backward_atom(value_type* in)
  {
    fftw_execute_dft_c2r(backward_plan,magle(in),in);
  }

  inline void transform_forward(value_type* in, int howmany=1)
  {
    value_type* offset=in;
    for(int i=0; i<howmany; ++i,offset+=fft_mem_offset) 
      transform_forward_atom(offset);
  }

  inline void transform_backward(value_type* in, int howmany=1)
  {
    value_type* offset=in;
    for(int i=0; i<howmany; ++i,offset+=fft_mem_offset) 
      transform_backward_atom(offset);
  }
};

/** specialization for multiple-fft with FFTW on complex
 */
template<typename T>
struct FFT_DFT<T,std::complex<T>,true,FFTW_ENG> 
: public fft_base, public fftw_base
{
  typedef std::complex<T> value_type;

  FFT_DFT(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    init(dims,in,howmany,uflag);
  }

  void init( int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    if(dims != fft_size)
    {
      fft_base::resize(dims,howmany);
#pragma omp critical
      {
        fftw_base::clean();
        fftw_iodim data_dims;
        fftw_iodim howmany_dims;
        data_dims.n=dims; data_dims.is=1; data_dims.os=1;
        howmany_dims.n=howmany; howmany_dims.is=dims; howmany_dims.os=dims;
        forward_plan = fftw_plan_guru_dft(1,&data_dims, 1,&howmany_dims
            , magle(in), magle(in),FFTW_FORWARD,uflag);
        backward_plan = fftw_plan_guru_dft(1,&data_dims, 1,&howmany_dims
            , magle(in), magle(in),FFTW_BACKWARD,uflag);
      }
    }
  }

  inline void transform_forward(value_type* in, int howmany=1)
  {
    fftw_execute_dft(forward_plan,magle(in),magle(in));
  }

  inline void transform_backward(value_type* in, int howmany=1)
  {
    fftw_execute_dft(backward_plan,magle(in),magle(in));
  }
};

/** specialization for multiple-fft with FFTW on double
 */
template<typename T>
struct FFT_DFT<T,T,true,FFTW_ENG> 
: public fft_base, public fftw_base
{
  typedef T value_type;

  FFT_DFT(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    init(dims,in,howmany,uflag);
  }

  void init(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    if(dims != fft_size)
    {
      fft_base::resize(dims,howmany);
#pragma omp critical
      {
        fftw_base::clean();

        fftw_iodim data_dims;
        data_dims.n=dims; data_dims.is=1; data_dims.os=1;

        fftw_iodim howmany_dims;
        howmany_dims.n=howmany; 
        howmany_dims.is=dims+2; howmany_dims.os=dims/2+1;
        forward_plan = fftw_plan_guru_dft_r2c(1,&data_dims, 1,&howmany_dims
            , in, magle(in),uflag);

        howmany_dims.is=dims/2+1; howmany_dims.os=dims+2;
        backward_plan = fftw_plan_guru_dft_c2r(1,&data_dims, 1,&howmany_dims
            , magle(in), in,uflag);
      }
    }
  }

  inline void transform_forward(value_type* in, int howmany=1)
  {
    fftw_execute_dft_r2c(forward_plan,in,magle(in));
  }

  inline void transform_backward(value_type* in, int howmany=1)
  {
    fftw_execute_dft_c2r(backward_plan,magle(in),in);
  }
};

//}
#endif
